This file recreates the statistics reported in Iosad, Pavel. 2016. The phonologization of redundancy: Length and quality in Welsh vowels. MS., University of Edinburgh. Submitted to Phonology.
This file is an R notebook (for more on the format see here). The code below reproduces the analyses, tables and figures in the table; you can use the notebook capabilities of RStudio to examine the data in detail, or you can download the RMarkdown source (under the ‘Code’ button in the top right corner).
The data is available as a separate R package from GitHub. To install this you can use the following:
library(devtools)
devtools::install_github('anghyflawn/llafaR')We now load the necessary packages. First we load the data itself, and make the data frame available as vowels.
library(llafaR)
data(vowels)These packages are needed to manipulate the data in various ways.
library(plyr)
library(dplyr)
library(car)
library(knitr)
library(tidyr)For plotting, we use ggplot2, and the package ggalt, which allows us to plot binned kernel density estimates.
library(ggplot2)
library(ggalt)Package mgcv provides generalized additive modelling capability.
library(mgcv)The following two packages are useful for reporting regression results, including the results of mgcv::gam.
library(texreg)
library(visreg)This package allows us to calculate the second-order AIC (AICc) of a model, which is recommended over standard AIC for our sample sizes.
library(MuMIn)The following are simply some customizations to simplify the visualization of data using ggplot2
theme_set(theme_bw())
vowel_lengths <- c('Short' = 'Short vowels', 'Long' = 'Long vowels')
vowel_cat_labeller <- as_labeller(function (x) paste0('//', x, '//'))The following function takes a list of models and produces a list of texreg objects suitable for presentation in table form; for conciseness, we exclude some goodness-of-fit statistics, but also add the AICc value to the table.
add_aicc_model_list <- function(model_list) {
llply(model_list, function (x) {
tr <- texreg::extract(x,
include.dispersion=F,
include.deviance=F,
include.dev.expl=F,
include.gcv=F,
include.nsmooth=F)
tr@gof <- c(AICc(x), tr@gof)
tr@gof.names <- c('AIC$_{c}$', tr@gof.names)
tr@gof.decimal <- c(TRUE, tr@gof.decimal)
return(tr)
})
}We are now ready to visualize the data for speaker Sp1. For the visualization alone, we do not need the data for [ə], as it does not enter the tenseness contrast. We also change the labels for the length categories
sp1.data <- vowels %>%
transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long'))) %>%
filter(speaker == 'Sp1', v1 != '\\sw')Now we calculate the mean duration and standard deviation of each or this speaker.
sp1.stats <- sp1.data %>%
group_by(v1, v1.is.long) %>%
summarise(v1h.dur.mean = mean(v1h.dur, na.rm=T),
v1h.dur.sd = sd(v1h.dur, na.rm = T),
v1h.dur.upr = v1h.dur.mean + v1h.dur.sd,
v1h.dur.lwr = v1h.dur.mean - v1h.dur.sd)We run a simple ANOVA to check if those duration figures are sufficiently different
sp1.dur.anova <- aov(v1h.dur ~ v1.is.long + v1, sp1.data)
summary(sp1.dur.anova) Df Sum Sq Mean Sq F value Pr(>F)
v1.is.long 1 0.3332 0.3332 554.36 < 2e-16 ***
v1 3 0.0428 0.0143 23.74 1.04e-13 ***
Residuals 278 0.1671 0.0006
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We now plot the distribution of duration by vowel quality and length category for Sp1
ggplot(sp1.data, aes(v1.is.long, v1h.dur * 1000))+
geom_violin()+
facet_wrap(~v1, labeller=vowel_cat_labeller)+
xlab('Vowel length')+
ylab('Duration including preaspiration, msec')+
geom_pointrange(data=sp1.stats,
aes(v1.is.long, v1h.dur.mean * 1000, ymin=v1h.dur.lwr * 1000, ymax=v1h.dur.upr * 1000),
colour='blue')This plot shows the vowel quality. Here we use colours, unlike the printed version.
ggplot(sp1.data, aes(v1.f2, v1.f1, colour=v1))+
geom_bkde2d()+
scale_x_reverse()+
scale_y_reverse()+
labs(x='F2', y='F1')+
facet_wrap(~v1.is.long,
labeller=as_labeller(vowel_lengths))+
scale_colour_brewer(palette='Set1',
name="Vowel",
labels=vowel_cat_labeller)We are now ready to do some stats. We do need the [ə] data for this purpose so we create a new data frame. We do recode the length categories for ease of further plotting.
sp1.model.data <- vowels %>%
filter(speaker == 'Sp1') %>%
mutate(v1.is.long = factor(recode(.$v1.is.long,
"T='Long';F='Short'"),
c('Short', 'Long')))Now we fit a range of GAMs. We start with almost the simplest possible model. It predicts the F1 of the stressed vowel from vowel quality and vowel duration; we also include a smooth of the F2.
sp1.fit.simple <- gam(v1.f1 ~ v1 + s(v1h.dur, k=5) + s(v1.f2, k=10),
data=sp1.model.data)
gam.check(sp1.fit.simple)
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 14 iterations.
The RMS GCV score gradient at convergence was 0.003052367 .
The Hessian was positive definite.
Model rank = 18 / 18
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(v1h.dur) 4.000 3.567 0.994 0.39
s(v1.f2) 9.000 8.812 1.052 0.80
Taking log2 of the response variable improves the residuals, so we stick to that
sp1.fit0 <- gam(log2(v1.f1) ~ v1 + s(v1h.dur, k=5) + s(v1.f2, k=10),
data=sp1.model.data)
gam.check(sp1.fit0)
Method: GCV Optimizer: magic
Smoothing parameter selection converged after 15 iterations.
The RMS GCV score gradient at convergence was 4.987601e-07 .
The Hessian was positive definite.
Model rank = 18 / 18
Basis dimension (k) checking results. Low p-value (k-index<1) may
indicate that k is too low, especially if edf is close to k'.
k' edf k-index p-value
s(v1h.dur) 4.00 3.53 1.03 0.77
s(v1.f2) 9.00 8.84 1.05 0.80
Now we fit a range of other models. First, we add the phonological length category.
sp1.fit <- update(sp1.fit0, . ~ . + v1.is.long)Then we add an interaction between vowel quality and phonological length
sp1.fit.inter <- update(sp1.fit, . ~ . + v1:v1.is.long)This model is the same as the previous one, but excludes the duration
sp1.fit.inter.nodur <- update(sp1.fit.inter, . ~ . - s(v1h.dur, k=5))Finally, we take the biggest model so far (the one with an interaction between quality and phonological length) and add a random intercept by word. We should do this anyway, because both vowel and length are nested within word. We model the random effects as a smooth component, because our random effect structure is very simple. (See here.).
sp1.fit.inter.random <- update(sp1.fit.inter, . ~ . + s(word, bs='re'))Here, we add the random effect to the model that has the vowel quality and length separately, but lacks the interaction.
sp1.fit.inter.random.nodur <- update(sp1.fit.inter.nodur, . ~ . + s(word, bs='re'))Finally, we add the random effects to the model without phonological length.
sp1.fit.nolength.random <- update(sp1.fit0, . ~ . + s(word, bs='re'))A quick look at \(AIC_c\) for all the models…
sp1_aicc_table <- AICc(sp1.fit0, sp1.fit, sp1.fit.inter, sp1.fit.inter.nodur,
sp1.fit.inter.random, sp1.fit.inter.random.nodur,
sp1.fit.nolength.random)
kable(sp1_aicc_table)| df | AICc | |
|---|---|---|
| sp1.fit0 | 18.37431 | -248.6609 |
| sp1.fit | 18.84652 | -270.0628 |
| sp1.fit.inter | 21.87908 | -322.7420 |
| sp1.fit.inter.nodur | 18.75322 | -302.2194 |
| sp1.fit.inter.random | 80.77597 | -337.9799 |
| sp1.fit.inter.random.nodur | 82.21295 | -327.9666 |
| sp1.fit.nolength.random | 94.43361 | -310.4034 |
Clearly, despite the much larger degrees of freedom, the models with the random effects are better than the ones without. We can now take a closer look at them. We look at the models with random effects, and just in case at the model with quality-length interaction but without the random effects. We need to run the argument to htmlreg through our add_aicc_model_list function in order to show \(AIC_c\). We also omit the effect for ‘long [ə]’, because that’s not actually found in the data. To establish significance, we can calculate 95% confidence intervals (if the CI excludes zero, the effect is significant). We could have used a log likelihood test (using R’s anova.gam) for the nested models as well.
htmlreg(add_aicc_model_list(list(sp1.fit.inter.random.nodur,
sp1.fit.nolength.random, sp1.fit.inter.random, sp1.fit.inter)),
caption="Models for $\\log_{2}$(F1) for speaker Sp1",
custom.model.names=c('No duration effect', 'No length effect', 'Both duration and length',
'No random effect'),
custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//",
"Long vowel", "longschwa", "Long //e//", "Long //o//", "Long //u//",
"F2 smooth" , "Random intercept for word", "Duration smooth"),
omit.coef="longschwa",
ci.force=T,
doctype=F,
bold=0.05,
indentation="")| No duration effect | No length effect | Both duration and length | No random effect | ||
|---|---|---|---|---|---|
| Intercept | 8.55* | 8.45* | 8.51* | 8.54* | |
| [8.41; 8.70] | [8.30; 8.59] | [8.37; 8.65] | [8.43; 8.66] | ||
| //ə// | 0.54* | 0.43* | 0.51* | 0.51* | |
| [0.33; 0.75] | [0.21; 0.64] | [0.31; 0.71] | [0.34; 0.67] | ||
| //e// | 0.76* | 0.70* | 0.77* | 0.77* | |
| [0.64; 0.88] | [0.61; 0.79] | [0.65; 0.89] | [0.68; 0.86] | ||
| //o// | 0.88* | 0.79* | 0.89* | 0.81* | |
| [0.61; 1.14] | [0.51; 1.07] | [0.63; 1.14] | [0.60; 1.03] | ||
| //u// | 0.12 | 0.20 | 0.10 | 0.02 | |
| [-0.16; 0.40] | [-0.09; 0.48] | [-0.17; 0.36] | [-0.20; 0.24] | ||
| Long vowel | -0.29* | -0.19* | -0.21* | ||
| [-0.46; -0.12] | [-0.36; -0.02] | [-0.35; -0.07] | |||
| Long //e// | -0.14 | -0.12 | -0.13* | ||
| [-0.30; 0.02] | [-0.27; 0.04] | [-0.25; -0.01] | |||
| Long //o// | -0.25* | -0.23* | -0.19* | ||
| [-0.47; -0.03] | [-0.44; -0.01] | [-0.37; -0.02] | |||
| Long //u// | 0.23 | 0.24* | 0.30* | ||
| [-0.02; 0.47] | [0.01; 0.47] | [0.11; 0.49] | |||
| F2 smooth | 8.77 | 8.68 | 8.76 | 8.71 | |
| [-8.72; 26.25] | [-8.72; 26.09] | [-8.72; 26.24] | [-8.88; 26.31] | ||
| Random intercept for word | 63.45 | 76.51 | 58.88 | ||
| [-146.27; 273.16] | [-141.05; 294.06] | [-150.84; 268.59] | |||
| Duration smooth | 3.24 | 3.14 | 3.17 | ||
| [-3.84; 10.33] | [-3.79; 10.07] | [-4.00; 10.34] | |||
| AIC\(_{c}\) | -327.97 | -310.40 | -337.98 | -322.74 | |
| AIC | -380.23 | -382.63 | -388.17 | -325.85 | |
| BIC | -64.24 | -19.67 | -77.70 | -241.76 | |
| Log Likelihood | 272.33 | 285.75 | 274.86 | 184.80 | |
| R2 | 0.95 | 0.95 | 0.95 | 0.93 | |
| Num. obs. | 345 | 345 | 345 | 345 | |
| * 0 outside the confidence interval | |||||
The biggest model is the most successful one.
Although it’s useful to know whether a particular effect is significant, interpreting them can be tricky, especially when they are non-linear or when there are interactions. It might be useful to visualize them. We use the package visreg for that. We create two visreg objects, for duration effect and for the effects of phonological length (which differ by vowel), and then visualize them using ggplot2 (visreg can do its own plotting, but we want consistency throughout the paper).
sp1.dur.effects <- visreg(sp1.fit.inter.random, xvar='v1h.dur', plot=F)
sp1.length.effects <- visreg(sp1.fit.inter.random, xvar='v1.is.long', by='v1', plot=F)First we plot the effect of duration. Since the original model has the response variable on the log scale, we need to exponentiate the fitted values back again. We also reverse the y-scale, because this is F1, so it makes sense if we want to visualize as the effect on vowel height.
ggplot(sp1.dur.effects$fit, aes(v1h.dur*1000, 2^visregFit))+
geom_line(colour='blue')+
geom_ribbon(aes(ymin=2^visregLwr, ymax=2^visregUpr), alpha=.2)+
geom_point(data=sp1.dur.effects$res, aes(v1h.dur*1000, 2^visregRes))+
xlab('Duration including preaspiration, msec')+
ylab('Estimated F1 (reversed scale)')+
scale_y_reverse()And now the phonological length effects. We don’t want the schwa here, because it’s never long.
ggplot(filter(sp1.length.effects$fit, v1 != '\\sw'), aes(v1.is.long, 2^visregFit))+
geom_jitter(data=filter(sp1.length.effects$res, v1 != '\\sw'),
aes(v1.is.long, 2^visregRes), size=1, alpha=.4)+
geom_crossbar(aes(ymin=2^visregLwr, ymax=2^visregUpr), colour='blue')+
facet_wrap(~v1, labeller=vowel_cat_labeller)+
xlab('Vowel length')+
ylab('Estimated F1 (reversed scale)')+
scale_y_reverse()First we filter out the speakers who do not show the south-western pattern.
sw.data <- filter(vowels, !(speaker %in% c('Sp1', 'Sp8')))We then prepare the data for visualization, much as with Sp1.
sw.data.vis <- sw.data %>%
transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long')))
sw.data.noschwa <- filter(sw.data.vis, v1 != '\\sw')
sw.stats <- sw.data.noschwa %>%
group_by(v1, v1.is.long) %>%
summarise(v1h.dur.mean = mean(v1h.dur.norm, na.rm=T),
v1h.dur.sd = sd(v1h.dur.norm, na.rm = T),
v1h.dur.upr = v1h.dur.mean + v1h.dur.sd,
v1h.dur.lwr = v1h.dur.mean - v1h.dur.sd)Again, an ANOVA to check that the difference in duration between short and long vowels is real
sw.dur.anova <- aov(v1h.dur.norm ~ v1.is.long + v1, sw.data.noschwa)
summary(sw.dur.anova) Df Sum Sq Mean Sq F value Pr(>F)
v1.is.long 1 242.1 242.09 744.9 <2e-16 ***
v1 3 119.0 39.67 122.0 <2e-16 ***
Residuals 1699 552.2 0.33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
First we plot the (normalized) duration by vowel quality and phonological length.
ggplot(sw.data.noschwa, aes(v1.is.long, v1h.dur.norm))+
geom_violin()+
facet_wrap(~v1, labeller=vowel_cat_labeller)+
xlab('Vowel length')+
ylab('Normalized duration including preaspiration')+
geom_pointrange(data=sw.stats, aes(v1.is.long, v1h.dur.mean, ymin=v1h.dur.lwr, ymax=v1h.dur.upr),
colour='blue')Then we plot the distribution of vowel qualities. We’re not interested in schwa, so that’s excluded here.
ggplot(sw.data.noschwa, aes(v1.f2.norm, v1.f1.norm, colour=v1))+
geom_bkde2d()+
scale_x_reverse()+
scale_y_reverse()+
xlab('Normalized F2')+
ylab('Normalized F1')+
facet_wrap(~v1.is.long, labeller=as_labeller(vowel_lengths))+
scale_colour_brewer(palette='Set1',
name='Vowel quality',
labels=vowel_cat_labeller)First we recode the word ffenestr from speakers other than Sp4, which is produced with [ɛː] instead of the expected [eː].
sw.data$v2.is.high[sw.data$speaker != 'Sp4' & sw.data$word == 'ffenestr'] <- TThe baseline model does not include any influence from the post-tonic vowel. It does include duration, quality and length (with interaction), and speaker and word random intercepts.
sw.fit.nohigh2 <- gam(v1.f1.norm ~ v1*v1.is.long + s(v1h.dur.norm, k=5)
+ s(v1.f2.norm, k=5) + s(speaker, bs='re') + s(word, bs='re'),
data=sw.data)Then we include a model with a single effect of post-tonic vowel height
sw.fit.high2.nointer <- gam(v1.f1.norm ~ v1*v1.is.long
+ v2.is.high + s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5)
+ s(speaker, bs='re') + s(word, bs='re'),
data=sw.data)The next model includes a three-way interaction between the quality and length of the stressed vowel and the height of the post-tonic vowel (this should allow to capture the described effect, where height only influences mid, long vowels).
sw.fit <- gam(v1.f1.norm ~ v1*v1.is.long*v2.is.high
+ s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5)
+ s(speaker, bs='re') + s(word, bs='re'),
data=sw.data)The model with the three-way interaction appears to be the best one.
aicc_list_sw <- AICc(sw.fit.nohigh2, sw.fit.high2.nointer, sw.fit)
kable(aicc_list_sw)| df | AICc | |
|---|---|---|
| sw.fit.nohigh2 | 117.6956 | 2113.372 |
| sw.fit.high2.nointer | 117.3676 | 2105.939 |
| sw.fit | 105.9340 | 2085.733 |
A more detailed look at the models shows that the model with the three-way interaction is the best of the three, and only shows significant effects of post-tonic vowel height if the stressed vowels are mid and long. This is in line with the descriptions.
htmlreg(add_aicc_model_list(list(sw.fit.nohigh2, sw.fit.high2.nointer, sw.fit)),
caption="Models for normalized F1, south-western speakers",
custom.model.names=c("No height effect", "No interaction", "Model with interaction"),
custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//", "Long vowel", "Long schwa", "Long /e/", "Long /o/", "Long /u/", "Duration smooth", "F2 smooth", "Speaker (random)", "Word (random)", "High post-tonic vowel", "schwa before high", "//e// before high", "//o// before high", "//u// before high", "Long vowel before high", "Long schwa before high", "Long //e// before high", "Long //o// before high", "Long //u// before high"),
omit.coef="schwa",
ci.force=T,
doctype=F,
indentation='',
bold=0.05)| No height effect | No interaction | Model with interaction | ||
|---|---|---|---|---|
| Intercept | -1.01* | -1.06* | -1.00* | |
| [-1.24; -0.77] | [-1.29; -0.83] | [-1.18; -0.83] | ||
| //ə// | 0.70* | 0.65* | 0.79* | |
| [0.44; 0.97] | [0.38; 0.91] | [0.58; 1.00] | ||
| //e// | 1.54* | 1.43* | 1.58* | |
| [1.28; 1.81] | [1.16; 1.70] | [1.35; 1.82] | ||
| //o// | 1.59* | 1.51* | 1.54* | |
| [1.27; 1.91] | [1.19; 1.83] | [1.26; 1.81] | ||
| //u// | 0.26 | 0.14 | 0.29 | |
| [-0.09; 0.60] | [-0.20; 0.49] | [-0.04; 0.62] | ||
| Long vowel | -0.22 | -0.28* | -0.25* | |
| [-0.49; 0.05] | [-0.56; -0.01] | [-0.47; -0.04] | ||
| Long /e/ | -0.26 | -0.17 | -0.83* | |
| [-0.62; 0.10] | [-0.52; 0.19] | [-1.15; -0.52] | ||
| Long /o/ | 0.00 | 0.06 | -0.38* | |
| [-0.35; 0.36] | [-0.29; 0.42] | [-0.68; -0.08] | ||
| Long /u/ | 0.34 | 0.34 | 0.35 | |
| [-0.09; 0.76] | [-0.09; 0.76] | [-0.16; 0.85] | ||
| Duration smooth | 1.86 | 1.53 | 2.23 | |
| [-2.70; 6.42] | [-2.17; 5.24] | [-3.16; 7.61] | ||
| F2 smooth | 3.34 | 3.48 | 3.79 | |
| [-4.05; 10.73] | [-4.06; 11.01] | [-3.97; 11.56] | ||
| Speaker (random) | 4.41 | 4.43 | 4.35 | |
| [-5.39; 14.21] | [-5.37; 14.23] | [-5.45; 14.15] | ||
| Word (random) | 98.08 | 96.93 | 76.56 | |
| [-117.51; 313.68] | [-118.67; 312.52] | [-123.35; 276.48] | ||
| High post-tonic vowel | 0.27* | 0.05 | ||
| [0.15; 0.38] | [-0.27; 0.36] | |||
| //e// before high | -0.08 | |||
| [-0.46; 0.29] | ||||
| //o// before high | 0.02 | |||
| [-0.35; 0.38] | ||||
| //u// before high | -0.18 | |||
| [-0.60; 0.24] | ||||
| Long vowel before high | 0.03 | |||
| [-0.35; 0.41] | ||||
| Long //e// before high | 1.06* | |||
| [0.58; 1.54] | ||||
| Long //o// before high | 0.82* | |||
| [0.35; 1.30] | ||||
| Long //u// before high | 0.05 | |||
| [-0.59; 0.69] | ||||
| AIC\(_{c}\) | 2113.37 | 2105.94 | 2085.73 | |
| AIC | 2098.96 | 2091.61 | 2074.11 | |
| BIC | 2761.47 | 2752.27 | 2670.42 | |
| Log Likelihood | -931.78 | -928.44 | -931.12 | |
| R2 | 0.79 | 0.79 | 0.79 | |
| Num. obs. | 2057 | 2057 | 2057 | |
| * 0 outside the confidence interval | ||||
Finally, we fit two models that attempt to derive the influence of the post-tonic vowel not from its phonological height but from its F1 or duration.
sw.fit.v2f1 <- gam(v1.f1.norm ~ v1*v1.is.long
+ s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5) +s(v2.f1.norm, k=5)
+ s(speaker, bs='re') + s(word, bs='re'),
data=sw.data)
sw.fit.v2dur <- gam(v1.f1.norm ~ v1*v1.is.long
+ s(v1h.dur.norm, k=5) + s(v1.f2.norm, k=5) +s(v2.dur.norm, k=5)
+ s(speaker, bs='re') + s(word, bs='re'),
data=sw.data)The model with the phonological height is superior to the other ones:
aicc_list_sw_posttonic <- AICc(sw.fit, sw.fit.v2f1, sw.fit.v2dur)
kable(aicc_list_sw_posttonic)| df | AICc | |
|---|---|---|
| sw.fit | 105.9340 | 2085.733 |
| sw.fit.v2f1 | 120.6421 | 2117.256 |
| sw.fit.v2dur | 119.2208 | 2112.814 |
We prepare the data for this speaker. We also recode those words where all tokens of short [ɔ] are produced as [ʊ], but not the others. Finally, we prepare a data frame for visualization, without the [ə]
sp8.data <- filter(vowels, speaker == 'Sp8')
sp8.data$v1[sp8.data$word %in% c('modryb', 'popeth', 'posib', 'copi', 'honni', 'tonnau')] <- 'u'
sp8.data.noschwa <- filter(sp8.data, v1 != '\\sw')The prepare a plot showing the duration of vowels and post-tonic consonants for this speaker, we need to get the duration that’s stored in two different variables. We use tidyr to get that data out and then calculate the means for vowel and consonant duration for each vowel and context.
sp8.data %>%
transform(v1.is.long=c('Short', 'Long')[as.factor(v1.is.long)]) %>%
transform(v1.is.long=factor(v1.is.long, c('Short', 'Long'))) %>%
filter(v1 != '\\sw') %>%
select(v1, v1.is.long, v1h.dur, c1.dur) %>%
gather(key, value, v1h.dur:c1.dur) %>%
mutate(is.v=c('Consonant', 'Vowel')[as.factor(key == 'v1h.dur')]) %>%
transform(is.v=factor(is.v, c('Vowel', 'Consonant'))) -> sp8.plot.data
sp8.plot.data %>%
mutate(value=value*1000) %>%
group_by(is.v, v1, v1.is.long) %>%
summarize(mean.dur = mean(value, na.rm=T),
sd.dur = sd(value, na.rm=T),
dur.upr = mean.dur + sd.dur,
dur.lwr = mean.dur - sd.dur) -> sp8.plot.meansWe are now ready to plot. The differences in vowel duration between ‘short’ and ‘long’ contexts are not dramatic. However, importantly, the duration of the consonants also differs between these contexts, much as in other South Welsh varieties.
ggplot(sp8.plot.data, aes(v1.is.long, 1000*value, fill=is.v))+
geom_violin(position=position_dodge(width=1))+
facet_wrap(~v1, labeller=vowel_cat_labeller)+
xlab('Vowel length')+ylab('Duration (including preaspiration for vowels), msec')+
geom_pointrange(data=sp8.plot.means,
aes(v1.is.long, mean.dur, ymin=dur.lwr, ymax=dur.upr, colour=is.v),
position=position_dodge(width=1))+
scale_fill_brewer(palette='Set1', name='Type')+
scale_colour_brewer(name='Type')An ANOVA shows that there are differences in duration between short and long contexts, although not for all vowels. This result needs verification.
sp8.anova <- aov(v1h.dur ~ v1 + v1.is.long + v1:v1.is.long, sp8.data.noschwa)
summary(sp8.anova) Df Sum Sq Mean Sq F value Pr(>F)
v1 3 0.12187 0.04062 50.375 < 2e-16 ***
v1.is.long 1 0.01443 0.01443 17.891 3.19e-05 ***
v1:v1.is.long 3 0.01550 0.00517 6.408 0.000328 ***
Residuals 276 0.22256 0.00081
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We can also look at the vowel qualities for this speaker
sp8.data.noschwa %>%
transform(v1.is.long=factor(c('Short', 'Long')[as.factor(v1.is.long)], c('Short', 'Long'))) %>%
ggplot(., aes(v1.f2.norm, v1.f1.norm, colour=v1))+
geom_bkde2d()+
scale_x_reverse()+
scale_y_reverse()+
xlab('Normalized F2')+
ylab('Normalized F1')+
facet_wrap(~v1.is.long)+
scale_colour_brewer(palette='Set1',
name='Vowel quality',
labels=vowel_cat_labeller)We fit three models with F1 as the dependent variable. The biggest model includes both duration and a quality-length interaction for the stressed vowel, and we compare that with models that either exclude the interaction (but keep both length and quality), and one where duration does not have an effect.
sp8.fit.nointer <- gam(v1.f1 ~ v1 + v1.is.long + s(v1h.dur, k=5) +
s(v1.f2, k=5) + s(word, bs='re'),
data=sp8.data)
sp8.fit.nodur <- gam(v1.f1 ~ v1*v1.is.long + s(v1.f2, k=5) +
s(word, bs='re'),
data=sp8.data)
sp8.fit <- gam(v1.f1 ~ v1*v1.is.long + s(v1h.dur, k=5, by=v1) +
s(v1.f2, k=5) + s(word, bs='re'),
data=sp8.data)A quick \(AIC_c\) comparison shows the model without the interaction to be best (by a small margin). The model without duration is clearly inferior to both models that do include it.
sp8_aicc_list <- AICc(sp8.fit.nodur, sp8.fit.nointer, sp8.fit)
kable(sp8_aicc_list)| df | AICc | |
|---|---|---|
| sp8.fit.nodur | 30.62959 | 4107.098 |
| sp8.fit.nointer | 32.83353 | 4090.267 |
| sp8.fit | 37.44546 | 4093.566 |
A closer look at the three models
htmlreg(add_aicc_model_list(list(sp8.fit.nodur, sp8.fit.nointer, sp8.fit)),
caption="Models for F1, Sp8",
custom.model.names=c("No duration effect", "Duration effect",
"Model with interaction"),
custom.coef.names=c("Intercept", "//ə//", "//e//", "//o//", "//u//",
"Long vowel", "Long schwa", "Long /e/", "Long /o/",
"Long /u/", "F2 smooth", "Word (random)", "Duration smooth",
"Duration for //i//", "Duration for schwa",
"Duration for //e//", "Duration for //o//",
"Duration for //u//"),
omit.coef="schwa",
ci.force=T,
doctype=F,
indentation='',
bold=0.05)| No duration effect | Duration effect | Model with interaction | ||
|---|---|---|---|---|
| Intercept | 345.73* | 381.57* | 345.09* | |
| [286.15; 405.32] | [327.99; 435.16] | [271.91; 418.28] | ||
| //ə// | 99.06* | 88.54* | 155.20* | |
| [33.17; 164.95] | [30.06; 147.02] | [67.97; 242.42] | ||
| //e// | 294.38* | 250.44* | 277.96* | |
| [234.12; 354.65] | [199.97; 300.90] | [203.94; 351.98] | ||
| //o// | 197.05* | 114.78* | 191.61* | |
| [93.96; 300.13] | [19.17; 210.40] | [79.28; 303.94] | ||
| //u// | 47.06 | 31.88 | 59.40 | |
| [-50.98; 145.11] | [-57.61; 121.37] | [-48.50; 167.30] | ||
| Long vowel | 13.43 | -12.86 | 8.24 | |
| [-44.22; 71.08] | [-40.49; 14.77] | [-51.54; 68.01] | ||
| Long /e/ | -2.26 | -19.40 | ||
| [-77.03; 72.51] | [-96.29; 57.50] | |||
| Long /o/ | -63.91 | -74.77 | ||
| [-143.35; 15.53] | [-157.25; 7.72] | |||
| Long /u/ | 4.19 | 37.05 | ||
| [-81.50; 89.88] | [-51.42; 125.53] | |||
| F2 smooth | 2.77 | 2.76 | 2.81 | |
| [-3.63; 9.17] | [-3.60; 9.11] | [-3.66; 9.28] | ||
| Word (random) | 17.86 | 21.02 | 16.89 | |
| [-193.82; 229.54] | [-196.53; 238.58] | [-194.79; 228.56] | ||
| Duration smooth | 2.05 | |||
| [-2.86; 6.96] | ||||
| Duration for //i// | 1.49 | |||
| [-2.02; 4.99] | ||||
| Duration for //e// | 1.00 | |||
| [-0.96; 2.96] | ||||
| Duration for //o// | 1.40 | |||
| [-1.87; 4.68] | ||||
| Duration for //u// | 2.86 | |||
| [-3.59; 9.31] | ||||
| AIC\(_{c}\) | 4107.10 | 4090.27 | 4093.57 | |
| AIC | 4100.75 | 4082.94 | 4083.92 | |
| BIC | 4217.76 | 4208.36 | 4226.97 | |
| Log Likelihood | -2019.75 | -2008.64 | -2004.52 | |
| R2 | 0.67 | 0.69 | 0.69 | |
| Num. obs. | 337 | 337 | 337 | |
| * 0 outside the confidence interval | ||||
The effect of duration in the best model goes in the opposite direction to the effect of duration for speaker Sp1.
sp8.vis <- visreg(sp8.fit.nointer, xvar='v1h.dur', plot=F)
ggplot(sp8.vis$fit, aes(v1h.dur * 1000, visregFit))+
geom_line(colour='blue', size=1)+
geom_ribbon(aes(ymin=visregLwr, ymax=visregUpr),
alpha=.3)+
geom_point(data=sp8.vis$res,
aes(v1h.dur * 1000, visregRes),
size=1, alpha=.3)+
scale_y_reverse()+
xlab('Duration including preaspiration, msec')+
ylab('Estimated F1 (reversed scale)')We can plot the vowel qualities in final syllables depending on whether there is a coda. There is a difference in high vowels but not in the other ones.
vowels %>% transform(final.closed=c('Open', 'Closed')[as.factor(final.closed)]) %>%
filter(v2 != 'h', !(is.na(final.closed))) %>%
ggplot(aes(v2.f2.norm, v2.f1.norm, colour = v2))+
geom_bkde2d()+
scale_x_reverse()+
scale_y_reverse()+
xlab('Normalized F2')+
ylab('Normalized F1')+
facet_wrap(~final.closed)+
scale_colour_brewer(palette='Set1',
name='Post-tonic vowel',
labels=vowel_cat_labeller)We verify that difference by building a model that includes an interaction between the presence of a coda and the quality of the vowel, and then plotting the estimated effect
posttonic.fit <- gam(v2.f1.norm ~ v2*final.closed + s(v2.dur.norm, k=5) +
s(v2.f2.norm, k=10) + s(speaker, bs='re'),
data=vowels)
pt.vis <- visreg(posttonic.fit, 'final.closed', 'v2', plot=F)
ggplot(pt.vis$fit, aes(final.closed, visregFit))+
geom_jitter(data=pt.vis$res,
aes(final.closed, visregRes),
size=1, alpha=.1)+
geom_crossbar(aes(ymin=visregLwr, ymax=visregUpr),
colour='blue')+
facet_wrap(~v2,ncol=3, labeller=vowel_cat_labeller)+
xlab('Final syllable type')+
ylab('Estimated normalized F1')+
scale_x_discrete(labels=c('Open', 'Closed'))